SURFACE COMPARISON WITH MASS TRANSPORTATION 



Y. LIPMAN, I. DAUBECHIES 

Abstract. We use mass-transportation as a tool to compare surfaces (2-manifoIds). In particular, we 
determine the "similarity" of two given surfaces by solving a mass-transportation problem between their 
conformal densities. This mass transportation problem differs from the standard case in that we require the 
solution to be invariant under global Mobius transformations. 

Our approach provides a constructive way of defining a metric in the abstract space of simply-connected 
smooth surfaces with boundary (i.e. surfaces of disk-type); this metric can also be used to define meaningful 
intrinsic distances between pairs of "patches" in the two surfaces, which allows automatic alignment of the 
surfaces. We provide numerical experiments on "real-life" surfaces to demonstrate possible applications in 
natural sciences. 



1. INTRODUCTION 

Alignment of surfaces plays a role in a wide range of scientific disciplines. It is a standard problem in com- 
paring different scans of manufactured objects; various algorithms have been proposed for this purpose in 
the computer graphics literature. It is often also a crucial step in a variety of problems in medicine and 
biology; in these cases the surfaces tend to be more complex, and the alignment problem may be harder. For 
instance, neuroscientists studying brain function through functional Magnetic Resonance Imaging (fMRI) 
typically observe several people performing identical tasks, obtaining readings for the corresponding activity 
in the brain cortex of each subject. In a first approximation, the cortex can be viewed as a highly convoluted 
2-dimensional surface. Because different cortices are folded in very different ways, a synthesis of the observa- 
tions from different subjects must be based on appropriate mappings between pairs of brain cortex surfaces, 
which reduces to a family of surface alignment problems fHl 127] . In another example, paleontologists studying 
molar teeth of mammals rely on detailed comparisons of the geometrical features of the tooth surfaces to 
distinguish species or to determine similarities or differences in diet j2] . 

Mathematically, the problem of surface alignment can be described as follows: given two 2-surfaces M and 
A/", find a mapping f : M. ^ M that preserves, as best possible, "important properties" of the surfaces. 
The nature of the "important properties" depends on the problem at hand. In this paper, we concentrate 
on preserving the geometry, i.e., we would like the map / to preserve intrinsic distances, to the extent 
possible. In terms of the examples listed above, this is the criterion traditionally selected in the computer 
graphics literature; it also corresponds to the point of view of the paleontologists studying tooth surfaces. 
To align cortical surfaces, one typically uses the Talairach method IBj (which relies on geometrically defined 
landmarks and is thus geometric in nature as well), although alignment based on functional correspondences 
has been proposed more recently [27j . 

In this paper we propose a procedure to "geometrically" align surfaces, based on uniformization theory 
and optimal mass transportation. This approach is related to the computer graphics constructions in jl8) . 
which rely on the representation of isometrics between topologically equivalent simply-connected surfaces by 
Mobius transformations between their uniformization spaces, and which exploit that 1) the Mobius group 
has small dimensionality (e.g. 3 for disk-type surfaces and 6 for sphere-type) and 2) changing the metric in 
one piece of a surface has little infiuence on the uniformization of distant parts. These two observations lead, 
in IH], to fast and particularly effective algorithms to identify near- isometrics between differently deformed 
versions of a surface. In our present context, these same observations lead to a simple algorithm for surface 
alignment, reducing it to a linear programming problem. 
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We shall restrict ourselves to (sufRciently smooth) disk-type surfaces; we map them to metric densities defined 
on the hyperbolic disk, their canonical uniformization space. (Apart from simplifying the description of the 
surface, this also removes any effect of global translations and rotations on the description of each individual 
surface.) The alignment problem can then be studied in the framework of Kantorovich mass-transportation 
[13] between these metric densities, as follows. Mass-transportation seeks to minimize the "average distance" 
over which mass needs to be "moved" (in the most efficient such moving procedure) to transform one mass 
density /x into another, v. In our case the uniformizing metric density (or conformal factor) corresponding 
to an initial surface is not unique, but is defined only up to a Mobius transformation. Because a naive 
application of mass-transportation on the hyperbolic disk would not possess the requisite invariance under 
Mobius transformations, we generalize the mass-transportation framework, and replace the metric d(x, y) 
traditionally used in defining the "average displacement distance" by a metric that depends on ji and v, 
measuring the dissimilarity between the two metric densities on neighborhoods of x and y. Introducing 
neighborhoods also makes the definition less sensitive to noise in practical applications. The optimal way of 
transporting mass in this generalized framework, in which the orientation in space of the original surfaces is 
"factored away" , automatically defines a corresponding optimal way of aligning the surfaces. 

Our approach also allows us to define a new distance between surfaces. The average distance over which 
mass needs transporting (to transform one metric density into the other) quantifies the extent to which the 
two surfaces differ; we prove that it defines a distance metric between surfaces. 

Other distances between surfaces have been used recently for several applications |19j . A prominent mathe- 
matical approach to define distances between surfaces considers the surfaces as special cases of metric spaces, 
and uses then the Gromov-Hausdorff (GH) distance between metric spaces [5]. The GH distance between 
metric spaces X and Y is defined through examining all the isometric embedding of X and Y into (other) 
metric spaces; although this distance possesses many attractive mathematical properties, it is inherently hard 
computationally |20[ IT]. For instance, computing the GH distance is equivalent to a non-convex quadratic 
programming problem; solving this directly for correspondences is equivalent to integer quadratic assign- 
ment, and is thus NP-hard [5]. In addition, the non-convexity implies that the solution found in practice 
may be a local instead of a global minimum, and is therefore not guaranteed to give the correct answer for 
the GH distance. The distance metric between surfaces that we define in this paper does not have these 
shortcomings: because the computation of the distance between surfaces in our approach can be recast as a 
linear program, it can be implemented using efficient polynomial algorithms that are moreover guaranteed 
to converge to the correct solution. 

It should be noted that in p^, Memoli generalizes the GH distance of j20j by introducing a quadratic mass 
transportation scheme to be applied to metric spaces already equipped with a measure (mm spaces); he 
notes that the computation of this Gromov-Wasserstein distance for mm spaces is somewhat easier and more 
stable to implement than the original GH distance. In our approach we do not need to equip the surfaces 
we compare with a measure: after uniformization reduces the problem to comparing two disks, we naturally 
"inherit" two corresponding conformal factors that we interpret as measure densities, for which we then 
apply an approach similar to the one proposed in |19j . Another crucial aspect in which our work differs from 
[l9] is that, in contrast to the (continuous) quadratic programming method proposed in [19j to compute the 
Gromov-Wasserstein distance between mm spaces, our conformal approach leads to a convex (even linear) 
problem, solvable via a linear programming method. 

It is worth mentioning that optimal mass transportation has been used as well, in the engineering literature 
to define interesting metrics between images; in this context metric is often called the Wasserstein distance. 
The seminal work for this image analysis approach is the paper by Rubner et al. J26], in which images are 
viewed as discrete measures, and the distance is called appropriately the "Earth Mover's Distance" . 

Another related method is presented in the papers of Zeng et al. [311 132] > which also use the uniformization 
space to match surfaces. Our work differs from that of Zeng et al. in that they use prescribed feature 
points (defined either by the user or by extra texture information) to calculate an interpolating harmonic 
map between the uniformization spaces, and then define the final correspondence as a composition of the 
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uniformization maps and this harmonic interpolant. This procedure is highly dependent on the prescribed 
feature points, provided as extra data or obtained from non-geometric information. In contrast, our work 
does not use any prescribed feature points, or external data, and makes use of only the geometry of the 
surface; in particular we make use of the conformal structure itself to define deviation from (local) isometry. 

Our paper is organized as follows: in Section[2]we briefly recall some facts about uniformization and optimal 
mass transportation that we shall use, at the same time introducing our notation. Section [3] contains the 
main results of this paper, constructing the distance metric between disk-type surfaces, in several steps. 



Section 
Section 



discusses various issues that concern the numerical implementation of the framework we propose; 
illustrates our results with a few examples. 



2. Background and Notations 



As described in the introduction, our framework makes use of two mathematical theories: uniformization 
theory, to represent the surfaces as measures defined on a canonical domain, and optimal mass transportation, 
to align the measures. In this section we recall some of their basic properties, and we introduce our notations. 



2.1. Uniformization. By the celebrated uniformization theory for Riemann surfaces (see for example 
lllp. any simply-connected Riemann surface is conformally equivalent to one of three canonical domains: 
the sphere, the complex plane, or the unit disk. Since every 2-manifold surface M equipped with a smooth 
Riemannian metric g has an induced conformal structure and is thus a Riemann surface, uniformization 
applies to such surfaces. Therefore, every simply- connected surface with a Riemannian metric can be 
mapped conformally to one of the three canonical domains listed above. We shall consider surfaces A4 that 
are topologically equivalent to disks and that come equipped with a Riemannian metric tensor g (possibly 
inherited from the standard 3D metric if the surface is embedded in R'^). For each such there exists a 
conformal map (f) : M ^ where V = {z \ \z\ < 1} is the open unit disk. The map pushes g to a metric 
on T); denoting the coordinates mV hy z — + ix^, we can write this metric as 

'g ^ <t>*g ~ ]i{z) Sij dx^ ® dx-' , 

where Jl{z) > 0, Einstein summation convention is used, and the subscript * denotes the "push-forward" 
action. The function Jl can also be viewed as the density function of the measure vol^vi induced by the 
Riemann volume element: indeed, for (measurable) A d M, 

(2.1) vo1a^(A)== / Jl{z) dx^ A dx'^ . 

It will be convenient to use the hyperbolic metric on the unit disk (1 — \z\'^)^'^5ijdx^ (8) dx^ as a reference 
metric, rather than the standard Euclidean Sijdx^ (E) dx^; note that they are conformally equivalent (with 
conformal factor (1 — jzp)"^). Instead of the density /i(z), we shall therefore use the hyperbolic density 
function 

(2.2) ^^"{z)■.^{l~\zffJl{z), 

where the superscript H stands for hyperbolic. We shall often drop this superscript: unless otherwise stated 
H — , and h' = ly^ in what follows. The density function /i = fi^ satisfies 



vol^ {A) = 




n{z)dvolH{z) , 



where dvol/f (z) = (1 — |zp) ^ dx^ A dx^. 

The conformal mappings of V to itself are the disk-preserving Mobius transformations m € Mjj, a family 
with three real parameters, defined by 

(2.3) m(z) = e'^^^, a G P, ^ £ [0, 2^). 
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Since these Mobius transformations satisfy 

(2.4) (l-|m(z)n-2|^'(^)|2 = (l_|^|2)-2^ 

where m! stands for the derivatives of m, the pull-back of /i under a mapping m e Md takes on a particularly 
simple expression. Setting w — to(z), with w = + iy^, and g{w) — Jl{w)Sijdy'^ (g) dy^ — ^{w){l — 
\w\'^)~'^Sijdy'^ (g) dy^ , the definition 

{m*'g){z)ki dx'' (g) dx^ := n{w) (1 - |w|^)"^ Sij dy' (g) dy^ 



implies 



{m*g)Miz) dx'' (g) dx*^ ^ fi{m{z)){l ^ \m{z)\^) l^c '^^'^ ® '^^^ 

= m("i(z)) (1 - |m(z)n-2 |m'(z)p 4^ «) dx'^ 
= fJ.{miz)) (1 - 4£ dx'' (g) dx'^. 



In other words, {m*g){z)ki dx'' (g) da;^ takes on the simple form m* ii{z) (1 — |zp)^^ 5ki dx'^ (g) dx^, with 

(2.5) m* ^{z) — ^{m{z)). 

Likewise, the push-forward, under a disk Mobius transform m{z) ~ w, of the diagonal Riemannian metric 
defined by the density function /i = /i^, is again a diagonal metric, with (hyperbolic) density function 
m*/i(w) = (TO*/i)^ (w) given by 

(2.6) m^.fi{w) = fi{m~^{w)). 

It follows that checking whether or not two surfaces A4 and Af are isometric, or searching for (near-) 
isometrics between A4 and Af, is greatly simplified by considering the conformal mappings from A4, Af to V: 
once the (hyperbolic) density functions fi and v are known, it suffices to identify m e Mu such that i>{'m{z)) 
equals fi{z) (or "nearly" equals, in a sense to be made precise). This was exploited in [18] to construct fast 
algorithms to find corresponding points between two given surfaces. 

2.2. Optimal mass transportation. Optimal mass transportation was introduced by G. Monge [H], and 
L. Kantorovich [iT. It concerns the transformation of one mass distribution into another while minimizing a 
cost function that can be viewed as the amount of work required for the task. In the Kantorovich formulation, 
to which we shall stick in this paper, one considers two measure spaces X, Y, a probability measure on each, 
/X G P{X), V G P(y) (where P(X),P(y) are the respective probability measure spaces on X and Y\ and 
the space n(/i, j/) of probability measures on X x y with marginals /i and v (resp.), that is, for A C X, 
5 C y, 7r(A X F) = A*(^) and 7r(X x B) = v{B). The optimal mass transportation is the element of n(/i, v) 
that minimizes J^^y dix,y)dTT{x,y), where d{x,y) is a cost function. (In general, one should consider an 
infimum rather than a minimum; in our case, X and Y are compact, d{-, •) is continuous, and the infimum 
is achieved.) The corresponding minimum, 

(2.7) Tiil^,i^)= inf / dix,y)dw{x,y), 

is the optimal mass transportation distance between fi and v, with respect to the cost function d{x, y). 

Intuitively, one can interpret this as follows: imagine being confronted with a pile of sand on the one hand 
(/i), and a hole in the ground on the other hand (— i^), and assume that the volume of the sand pile equals 
exactly the volume of the hole (suitably normalized, fj,, v are probability measures). You wish to fill the hole 
with the sand from the pile (tt S n(/i, i^)), in a way that minimizes the amount of work (represented by 
/ d{x,y)dTT(x,y), where c?(-, •) can be thought of as a distance function). In the engineering literature, the 
distance T^{fi, v) is often called the "earth mover's distance" [26], a name that echoes this intuition. 

In what follows, we shall apply this framework to the density functions /i and v on the hyperbolic disk T) 
obtained by conformal mappings from two surfaces A^, A/", as described in the previous subsection. 
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The main obstacle to applying the Kantorovich transportation framework directly is that the density /i, char- 
acterizing the Riemannian metric on T) obtained by pushing forward the metric on Ai via the uniformizing 
map (/) : — >■ I?, is not uniquely defined: another uniformizing map (j)' : M ^ D may well produce a differ- 
ent fj! . Because the two representations are necessarily isometric ((/)~^ o maps M. isometrically to itself), 
we must have fj.'{m{z)) — ij,{z) for some m G Md- (In fact, m — cj)' o (p^^.) In a sense, the representation of 
(disk-type) surfaces Ai as measures over T) should be considered "modulo" the disk Mobius transformations. 

We thus need to address how to adapt the optimal transportation framework to factor out this Mobius 
transformation ambiguity. This is done by designing a special distance (or cost) functional d^j^{z,w) that 
depends on the conformal densities /i and u representing the two surfaces. (A fairly simple argument shows 
that a cost function that does not depend on ^ and u allows only trivial answers, such as d{z, w) =Q for all 
z, w.) As we shall see in the next section, this cost function will have an intuitive explanation: d^ ^{z, w) will 
measure how well an _R-sized neighborhood of z with density n can be matched isometrically to an i?-sized 
neighborhood of w with density v by means of a disk Mobius transformation. 



3. Optimal volume transportation for surfaces 



We want to measure distances between surfaces by using the Kantorovich transportation framework to 
measure the transportation between the metric densities on V obtained by uniformization applied to the 
surfaces. The main obstacle is that these metric densities are not uniquely defined; they are defined up to a 
Mobius transformation. In particular, if two densities ^ and v are related by = rn*/x (i.e. fi{z) — ;/(m(z))), 
where m G Mjj, then we want our putative distance between fi and to be zero, since they describe isometric 
surfaces, and could have been obtained by different uniformization maps of the same surface. A standard 
approach to obtain quantities that are invariant under the operation of some group (in our case, the disk 
Mobius transformations) is by minimizing over the possible group operations. For instance, we could set 



Distance(/x, i^) = inf inf / d(z,w) dTT{z,w) 



where n(/^, v) is the set of probability measures on I? x 2? with marginals fxvoln and vvoIh- In order for this 
to be computationally feasible, we would want the minimum to be achieved in some m, which would depend 
on /X and of course; let's denote this special minimizing m e AId by m^^^. This would mean 



Distance(/x, z^) = inf / d{z,w) dTr{z,w) 



(3.1) = inf / d(mp ,^(z), w) (i7r(z, ui) . 

If v were itself already equal to m'^/i, for some m' S M^, then we would expect the minimizing Mobius trans- 
formation to be TO^jy = m'; for tt supported on the diagonal d — {(z, z) ; z € V} d V xV, defined by tt{A) ~ 
Jj^ v{w) dvolniw), with A2 = {w; {w, w) G A}, one would then indeed have J^^v '^("^m,!^('^)' dTr{z, w) — 0, 



leading to Distance(/i, m'^/i) = 0. From (3.1) one sees that this amounts to using the same formula as for 



the standard Kantorovich approach with just one change: the cost function depends on /i and i'. 

We shall use a variant on this construction, retaining the principle of using cost functions d{-,-) in the 
integrand that depend on fi and v, without picking them necessarily of the form d{mfj_^^{z),w). In addition 
to introducing such a dependence, we also wish to incorporate some robustness into the evaluation of the 
distance between (or dissimilarity of) /z and h'. We shaU do this by using a cost function d^^,{z,w) that 
depends on a comparison of the behavior ^ and v on neighborhoods of z and w, mapped by m ranging over 
Mjj. The next subsection shows precisely how this is done. 

3.1. Construction of d-^^,{z,w). We construct djj^{z,'w) so that it indicates the extent to which a neigh- 
borhood of the point z in (2?, /i), the (conformal representation of the) first surface, is isometric with a neigh- 
borhood of the point w in {V, i/), the (conformal representation of the) second surface. We will need to define 
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two ingredients for this: the neighborhoods we wiU use, and how we shall characterize the (dis) similarity of 
two neighborhoods, equipped with different metrics. 

We start with the neighborhoods. 

For a fixed radius i? > 0, we define i^zo,R to be the hyperbolic geodesic disk of radius R centered at zq- 
The following gives an easy procedure to construct these disks. If zo = 0, then the hyperbolic geodesic disks 
centered at = are also "standard" (i.e. Euclidean) disks centered at 0: fio.fl = {z; \z\ < tb}, where 
rji = arctanh(r) = R. The hyperbolic disks around other centers are images of these central disks under 
Mobius transformations (= hyperbolic isometrics): setting m{z) = (z — zo)(l — zzq)^^, we have 



(3.2) 



zo,R 



simply rotates f2o,i?, around its 



If m', m" are two maps in Md that both map zq to 0, then m" o (m 
center, over some angle 9 determined by m' and to". From this observation one easily checks that (3.2 1 holds 
for any m € Mo that maps zq to 0. In fact, we have the following more general 

Lemma 3.1. For arbitrary z,w (zV and any R > 0, every disk Mobius transformation to G Mjj that maps 
z to w (i.e. w = m{z)) also maps flz.R to Qw,R- 



Next we define how to quantify the (dis)similarity of the pairs (yizo.R j /^) snd {^wo,R i Since (global) 
isometrics are given by the elements of the disk-preserving Mobius group Mu, we will test the extent to 
which the two patches are isometric by comparing {flwo,R. j ^ ) with all the images of {flzo^R , fJ- ) under Mobius 
transformations in M^ that take zq to wq. 

To carry out this comparison, we need a norm. Any metric gij{z)dx^ ® dx^ induces an inner product on 
the space of 2-covariant tensors, as follows: if a(z) = aij{z) dx^ (S) dx^ and b(z) — bij(z) dx^ (8) dx^ are two 
2-covariant tensors in our parameter space V, then their inner product is defined by 

(3.3) (a(z),b(z)) =a,,(z)6,,(z)5*'=(z)g^-^(z) ; 

as always, this inner product defines a norm, ||a||^ ~ aij{z) akg{z) g'^^{z) g^^{z). 

Now, let us apply this to the computation of the norm of the difference between the local metric on one 
surface, gij{z) = /x(z)(l — \z\'^)~^Sij, and hij{w) — i^{w){l — \w\'^)~^Sij, the pull-back metric from the other 
surface by a Mobius transformation to. Using (3.3 1,(2.5 1, and writing S for the tensor with entries Sij, we 
have: 



m(z) - i^(TO(z)))'(l - |zn-4 S,, Ski g''{z) g^\z) 



\\^,{z)il-\z\')-'^-,.{miz)){l-\zn 
1 



= 1- 



v(m{zY) 
H{z) 



We are now ready to define the distance function d^^{z^w): 



(3.4) 



inf 

meMn ,m{zo)=wo Jq 



I ^(z) - {m*iy){z) I dvolniz), 



where dvol/f (z) = (1 — |zp) ^ dx A dy is the volume form for the hyperbolic disk. The integral in (3.4) can 
also be written in the following form, which makes its invariance more readily apparent: 



(3.5) 



1 



v(m(z)) 



^i{z) 



dvoljiiiz) = 



11^ - m*i'\\z dvol7n(z). 



where dvol7vi(z) ~ Ai(z)(l — |zp) ^ dx^ A dx^ ~ \/\gij \ dx^ A dx^ is the volume form of the first surface Ai. 



The next Lemma shows that although the integration in (3.5) is carried out w.r.t. the volume of the first 
surface, this measure of distance is nevertheless symmetric: 
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Lemma 3.2. If m £ Mjj maps zq to wq, m(zo) — wq, then 



fi(z) — m*i'{z) dvoljj{z) 



Proof. By the pull-back formula (2.5), we have 



— m*i'{z) dvol/f (z) = 



fi{z) — h'{m{z)) dvol/f (z) 



Performing the change of coordinates z = m ^{w) in the integral on the right hand side, we obtain 



/m(O.Q__R.) 

where we have used that is a n isometry and therefore preserves the volume eleme nt d volg(w) = (1 — 
dy^ A dy'^. By Lemma 3.1 ■m{^zo,R) = ^wo,R ; using the push-forward formula (2.6| then allows to 
conclude. □ 

Note that our point of view in defining our "distance" between z and w differs from the classical point of 
view in mass transportation: Traditionally, d{z,w) is some sort of physical distance between the points z 
and w; in our case dj}^{z,'w) measures the dissimilarity of (neighborhoods of) z and w. 

The next Theorem lists some important properties of c?^^; its proof is given in Appendix A. 
Theorem 3.3. The distance function dj}^{z,w) satisfies the following properties 



(1) d^i,,,m'^Am^^izo),m2^{wo)) = d^^^izo,wo) 



Invariance under (well-defined) 
Mobius changes of coordinates 



(2) dj^'^^{zQ,wo) ^ d^^^{wo,zo) Symmetry 

(3) d^^^{zQ,wo) >0 Non-negativity 

(4) dj^MizojiJUo) = ^zo,R in [T^ilA ^-^d ^wo,R in {'C'tv) are isometric 

(5) d^j.^j,(m"i(zo),zo) = Reflexivity 

(6) 0^^1,^3(21, 23) < ^^1,^2(^1' ^2) + rf/?2,M3(^2, 23) Triangle inequality 



In addition, the function d^^ : T) xV — >■ K is continuous. To show this, we first look a little more closely at 
the family of disk Mobius transformations that map one pre-assigned point zq G V to another pre-assigned 
point Wq G 2?, over which one minimizes to define d^{zQ,'Wo). 

Definition 3.4. For any pair of points zq, wq E T), we denote by M d,zo,wo the set of Mobius transformations 
that map zq to wq. 

This family of Mobius transformations is completely characterized by the following lemma: 

Lemma 3.5. For any zo,wq e V, the set M£)^zo,wq constitutes a 1 -parameter family of disk Mobius trans- 
formations, parametrized continuously over (the unit circle). More precisely, every m G M]j_za,wa is of 
the form 

ro r\ / ^ --n, f \ za-woa l-z^WQa 

io.b) m[z)=T with a = a[ZQ,wo,a):= — ^ and r = t(2o, Wq, cr) := cr ^ — , 

1 — az 1 — ZgWQa 1 — ZqWoCt 

where a £ Si := {z E C ] |z| = l} can he chosen freely. 
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Proof. By (2.3 1, the disk Mobius transformations that map zq to aU have the form 



/ \ ill, ^ ~ ^0 ,1 . f 1,- -u • -1 / \ w + e^'^ZQ 

m,!, zrAz) = e ^ , the inverse ot which is m.,, , (w) = e ^ , 

"^'^"^ ' l + e-"/'z^u'' 

where e M can be set arbitrarily. It follows that the elements of M£, zq,wq are given by the family 
^7 Lo °"^i/'-^o: with 7 e M. Working this out, one finds that these combinations of Mobius transformations 



take the form (3.61, with a = e^i^-^\ □ 



We shall denote by mzo,wo,cr the special disk Mobius transformation defined by (3.6|. In view of our interest 
in d^\^, we also define the auxiliary function 

<^> -.VxVx Si — >C 

by $(zo, wojc) = /n(2„ /j) I ^^{'^) ~ ^(™2o,iuo,cr(-^)) I dvolniz). This function has the following continuity prop- 
erties, inherited from jj, and v: 

Lemma 3.6. 

• For each fixed {zq,wq), the function ^{zo,wq, •) is continuous on Si. 

• For each fixed a G Si, $(•, •, a) is continuous on T> x D. Moreover, the family ( $(•, •, a) J is equicon- 
tinuous. 

Proof. The proof of this Lemma is given in Appendix A. □ 



Note that since S^ is compact. Lemma 3.6 implies that the infimum in the definition of d^\, can be replaced 
by a minimum: 



d „(zo,t«o)= min / \ ^i{z) - i^{m{z)) {dvolniz) . 

m{zo)=wo Jn^g_R 

We have now all the building blocks to prove 

Theorem 3.7. If jj, and v are continuous from T) to R, then dj}^{z,w) is a continuous function onT) xT). 
Proof. Pick an arbitrary point {zq,wq) CzV xV, and pick e > arbitrarily small. 



By Lemma 3.6 there exists a, S > such that, for \zq — zo\ < S, \w[^ ~ wo\ < S , we have 

\^{zo,wo,cr)-<i>{z'Q,w'Q,a) \ <e, 
uniformly in a. Pick now arbitrary Zq, Wq so that \zo — Zq\, \wo — Wq\ < S. 

Let mzg,wo,a, resp. rnz'^.w'^^.a' be the minimizing Mobius transform in the definition of dj} ^{zo,wo), resp. 
d^,^{zo,WQ) ^ ^{zo,wo,cr) and d^[^{z'„,Wo) ^ ^{zo,W(),a') . 

It then follows that 
d^^^{zo,wo) = mm$(zo, wo,t) < $(zo,wo,cr') 

< ^{z'o,w'o,(7') + \<i>{zQ,WQ,a') - <i>{z'o,w'o,cr')\ = d^^^{z'o,w'o) + \^{zo,wo,a') - ^{z'o,w'o,a')\ 

< d^A^Q^^'a) + sup \^{zo,wa,uj) - ^(z^,,^^,^)! < d^_^(z[„ w^) + £ . 

Likewise d^,^{z'Q,w'„) < d^,,(zo, wq) + £, so that |d^„(zo, wq) - c?/?,t.(2o, w^o) | < £• □ 
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3.2. Incorporating d^^{z^'w) into the transportation framework. The next step in constructing the 
distance operator between surfaces is to incorporate the distance dj},^{z, w) defined in the previous subsection 
into the (generahzed) Kantorovich transportation model: 



(3.7) Ti^(M,^)= inf / d;^,,{z,w)dn{z,w). 

The main result is that this procedure (under some extra conditions) furnishes a metric between (disk-type) 
surfaces. 

Theorem 3.8. There exists tt* G n(/i, i^) such that 

d^ ^,{z,w)dTT* (zjw) = inf / d^ ^{z,w)dTT{z,w). 

Proof. This proof follows the same argument as in |30j . adapted here to our generaliz ed s etting. It uses the 

continuity of the distance function to derive the existence of a global minimum of (3.7). Let (Tr^ ) e 

1 I V /feeN 

n(/i, v) be a minimizer sequence of (3.7 1, for example by taking 



d^^^{z,w)dTTk{z,w) < T^ifi,!^) + ^■ 



Then this sequence of measures is tight, that is, for every e > 0, there exists a compact set C C V x V such 
that 7rfc(C) > 1 — e, for all fc e N. To see this, note that since V is separable and complete, the measures ^, 
V are tight measures (see [22] )• This means that for arbitrary e > 0, there exist compact sets A. i? C I? so 
that [i[A) > 1 - e/2 and v{B) > 1 - e/2. It then follows that, for all fc G N, 

^fc(A X B) = ^k{A xV)-nkiAx {V\B))> fi{A) - v{V \ B) = /i(A) - (1 - v{B)) >l-e. 



Since the set C = AxB(Z'Dx'D\s compact, this proves the claimed tightness of the family ( tt^ . By 

V / fcGN 

Prohorov's Theorem [22], a tight family of measures is sequentially weakly compact; in our case this means 
that ( TTfe ) has a weakly convergent subsequence ( tt^ ) ; by definition, its weak limit tt* satisfies, for 

V /feGN V "/riGN 

every bounded continuous function f onV xV, 

f{z,w)dTTk^{z,w)^ / f(z,w)dTT*{z,w). 
VxV JvxV 

Therefore, taking in particular the continuous fimction f{z,w) = dj^' ^{z,w), we obtain 
Tdit^i^)^ 1™ / f{z,w)d-Kk^{z,w) = / f{z,w)dTT*{z,w). 

"■^°°JVxV JvxV 

□ 



Under rather mild conditions, the "standard" Kantorovich transportation (2.7) on a metric spaces {X,d) 
defines a metric on the space of probability measures on X . We will prove that our generalization defines 
a distance metric as well. More precisely, we shall prove first that 

d«(x,AA) = r«(/i,^.) 

defines a semi-metric in the set of all disk-type surfaces. We shall restrict ourselves to surfaces that are 
suSiciently smooth to allow uniformization, so that they can be globally and conformally parameterized over 
the hyperbolic disk. Under some extra assumptions, we will prove that is a metric, in the sense that 
d^{A4, J\f) = implies that A4 and Af are isometric. 

For the semi-metric part we will again adapt a proof given in [30^ to our framework. In particular, we shall 
make use of the following "gluing lemma" : 
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Lemma 3.9. Let /Lti,/i2,M3 he three probability measures on V, and let 7ri2 G Il{fii, fj,2), tt23 G n(/i2,M3) 
be two transportation plans. Then there exist a probability measure tt on V x V x V that has 7ri2,7r23 as 
marginals, that is /^^^^ d7r(zi, Z2, 2^3) = dT:i2{zi, Z2), and J_^^^^ dT:{zi, Z2, z^) = (^7723(2:2,23). 

This lemma will be used in the proof of the following: 

Theorem 3.10. For two disk-type surfaces M — (V,^), J\f — (V,i^), let d^{M.,Af) he defined by 

Then defines a semi-metric on the space of disk-type surfaces. 

Proof. The symmetry of d^ implies symmetry for T^, by the following argument: 

inf / d^^,^{z,w)dT:{z,w) = irif / d^:{w,z)dTr{z,w) 

= inf / d^'{w,z)d'TT(w,z), where we have set tt{w, z) — tt{z,w) 
= Tf{v, fi) . ( use that tt G n(^, v) ^ t: e n(z^, ^)) 



The non-negativity of •) automatically implies T^{fi, v) > 0. 

Next we show that, for any Mobius transformation m, rj^(/i,m*/x) — 0. To see this, pick the transportation 
plan TT G n(/i,TO*/i) defined by 



f{z,w)dTT{z,w)— / f{z,m{z))fi{z)dvolH{z). 
VxV Jv 



On the one hand tt G n(/i, m^/i), since 



dT:{z,w) — / fi{z)dvolH{z), 
AxV J A 



o?7r(z, w)= / XB{w)dTT{z,w) 



and 



XBim{z))^{z)dvolH{z) = / XB{w)^J.4w)dvolHiw), 
V Jv 

where we used the change of variables w = m(z) in the last step. Furthermore, 7r(z, w) is concentrated on 
the graph of m, i.e. on {(z,m(z)) ; z £T>} C T> x V. Since dj}^^^^{z,m{z)) — for all z G P we obtain 
therefore Td{fi, m^fi) < J^^^ dj}m_^^{z, w)dTT{z, w) = 0. 

Finally, we prove the triangle inequality T^{^i,^3) < Tj*(^i,^2) + T)f (a^2, A^s) ■ To this end we follow the 
argument in the proof given in [30j (page 208). This is where we invoke the gluing Lemma stated above. 



We start by picking arbitrary transportation plans tti2 G n(/Lti,/X2) and 7r23 G n(/i2,/X3). By Lemma 3.9 
there exists a probability measure t: on T> x T> x T> with marginals 7ri2 and 7r23. Denote by tti^ its third 
marginal, that is 

d7r(zi, Z2, Z3) = d7ri3(zi, Z3). 

Z2ev 
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Then 



?d(Aii,M3)< / d^^^^^{zi,Z3)dT:i3{zi,Z3) = d^^^^^{zi, Z3)dn{zi, Z2, Z3) 

JtixV Jvx'Dx'D 



< 



VxVxV 



-/ d^^^^^{zi,Z2)dT:{zi,Z2,Z3) + d^^^^^{z2, Z3)dTr{zi, Z2, Z3) 

Jvx'Dx'D Jdx'Dx'D 

</ d^'^^^^{zi,Z2)dTri2{zi,Z2) + d^^,^^(z2,Z3)d7r23 (22,23), 

'Dx'D J'Dx'D 



where we used the triangle-inequahty for d^ ^ hsted in (Theorem 3.3). Since we can choose 7ri2 and 7r23 to 
achieve arbitrary close values to the infimum in eq. (3.7) the triangle inequality follows. □ 

To qualify as a metric rather than a semi-metric, (or T^) should be able to distinguish from each other 
any two surfaces (or measures) that are not "identical" , that is isometric. To prove that they can do so, 
we need an extra assumption: we shall require that the surfaces we consider have no self-isometries. More 
precisely, we require that each surface M that we consider satisfies the following definition: 

Definition 3.11. A surface M is said to be a singly g-nfittable (where e M, Q ^ 0) if, for all R > g, and 
all z e 2?, there is no other Mobius transformation rn other than the identity for which 

|/i(z) — /i(m(z))| dvol/f (z) = 0. 

R 

Remark 3.12. This definition can also be read as follows: Ai is singly g-nfittable if and only if, for all R> g, 
any two conformal factors fii and fi2 for Ai satisfy: 

(1) For all z ^ D there exists a unique minimum to the function w 0?^^ 112^^' 

(2) For all pairs {z,w) G I? x I? that achieve this minimum there exists a unique Mobius transformation 
for which the integral in (3.4) vanishes (with /ii in the role of fj,, and fi2 in that of v). 

Essentially, this definition requires that, from some sufficiently large (hyperbolic) scale onwards, there are 
no isometric pieces within {V, /i) (or (2?, v)). 

We start with a lemma, and then prove the main result of this subsection. 

Lemma 3.13. Let tt e n(/i, ;/) be such that J-^^j, ^^1/(2, w) d7r(z, w) = 0. Then, for all zo E T> and S > 0, 
there exists at least one point z e ^zo.s such that djj\^{z,w) = for some uo eT). 
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Proof. By contradiction: assume that there exists a disk ^za.s such that d^j^{z,w) > for all z G ^zq.s and 
all w e P. Since 



d7r(z,w) = / ii{z) dvolniz) > , 
a{zo,5)xv Jn{zo,s) 

the set n{zo,6) x I? contains some of the support of tt. It follows that 

d^^^{z,w)dTT{z,w) > , 

fl{zo,S)x'D 

which contradicts 

d^j^{z,w)dTi{z,w)< i d^,^{z,w)dTT{z,w)=0 
n{zo,5)xV ' Jvxv 



□ 



Theorem 3.14. Suppose that A4 and Af are two surfaces that are singly g-ufittahle. If <l^{M.,I\f) = 
for some R > g, then there exists a Mobius transformation m € Md that is a global isometry between 
M. = (fj/i) and N = {T>,v) (where /i and v are conformal factors of M and N , respectively). 

Proof. When d^{MM) = 0, there exists (see |30]) tt € n(^, v) such that 



VxV 



da,,.(.Z,w)dTT{z,w) = 0. 



Next, pick an arbitrary point zq & T> such that, for some wq € V, we have dj}j,{zo,wo) = 0. (The existence 



of such a pair is guaranteed by Lemma 3.13 ) This implies that there exists a unique Mobius transformation 
nio G Md that takes zq to wq and that satisfies v{mo{z)) — fi{z) for all z G ^zo,R- We define 

p* = sup{p ; ^^^^(zo, Wo) = 0}; 

clearly p* > R. The theorem will be proved if we show that p* = oo. We shall do this by contradiction, i.e. 
we assume p* < oo, and then derive a contradiction. 

So let's assume p* < oo. Consider flzo,p', the hyperbolic disk around zq of radius p* . (See Figure [l] for 
illustration.) Set e = (i? — g)/2, and consider the points on the hyperbolic circle C = d^lzo,p*-g-e- For 



every zi G C, consider the hyperbolic disk ^zi,e/2', by Lemma 3.13 there exists a point Z2 in this disk and a 
corresponding point W2 E T) such that d^ j^{z2,W2) — 0, i.e. such that 



\fi{z) - m'*iy{z)f dvolniz) = 

for some Mobius transformation m' that maps Z2 to W2', in particular, we have that 
(3.8) p.{z) — i>{m'{z)) for all z G ^Z2,R ■ 

The hyperbolic distance from Z2 to dVlz^.p' is at least g + e/2. It follows that the hyperbolic disk ^z'2,e+e/i 
is completely contained in fi^^ p*; since pl{z) = u{mo{z)) for all z G ri^^^p*, this must therefore hold, in 
particular, for all z G f^22,e+e/4- Since ^Z2.e+e/i C f^z2,-R' also have fi{z) — v{m'{z)) for all z G f^22,g+e/4, 



by (3.8). This implies v(w) = v{mQ o (m') ^(w)) for all w G f^u)2,e+e/4- Because M is singly g-nfittable, it 



follows that TTiQ o (m') ^ must be the identity, or mo = m'. Combining this with (3.8), we have thus shown 
that p{z) — v{mo{z)) for all z G ^Z2,R- 

Since the distance between Z2 and zi is at most e/2, we also have 

^Z2,R. ^ ^2i,_R-e/2 = ^zi,e+3e/2 ■ 

This implies that if we select such a point 22(^1) for each zi G C, then 51zf,.p*-g-£ U ( \Jz-^ec ^z2{zi).r) covers 
the open disk ^zo.p*+e/2- By our earlier argument, p,{z) = v{mo{z)) for all z in each of the ^z2(zi).r] since 
the same is true on ^zo,p* -g-e^ it follows that p,{z) = v(mQ{z)) for all z in i^zo,p*+e/2- This contradicts the 



Surface Comparison with Mass Transportation 



13 



definition of p* as the supremum of all radii for which this was true; it follows that our initial assumption, 
that p* is finite, cannot be true, completing the proof. □ 

For iV, p) to be singly g-nfittable, no two hyperbolic disks ^z,r, ^w,r (where w can equal z) can be isometric 
via a Mobius transformation m, if i? > g, except if m = Id. However, if z is close (in the Euclidean sense) 
to the boundary of P, the hyperbolic disk Q,z^r is very small in the Euclidean sense, and corresponds to 
a very small piece (near the boundary) of fA. This means that single gi-nfittability imposes restrictions in 
increasingly small scales near the boundary of A^; from a practical point of view, this is hard to check, and in 
many applications, the behavior of M close to its boundary is irrelevant. For this reason, we also formulate 
the following relaxation of the results above. 

Definition 3.15. A surface is said to be a singly ^-^fittable (where ^4 > 0) if there are no patches (i.e. 
open, path-connected sets) in M of area larger than A that are isometric, with respect to the metric on 7W. 

If a surface is singly yl-^fittable, then it is obviously also ^'-^fittable for all A' > A; the condition of 
being A-^fittable becomes more restrictive as A decreases. The following theorem states that two singly 
A-j^^fittable surfaces at zero d -distance from each other must necessarily be isometric, up to some small 
boundary layer. 

Theorem 3.16. Consider two surfaces M andAf, with corresponding conformal factors p and v on V, and 
suppose d^{A4,Af) — for some R > 0. Then the following holds: for arbitrarily large p > 0, there exist 
a Mobius transformation m G Mjj and a value A > such that if A4 and Af are singly A-^fittable then 
p{m{z)) = ^{z), for all z € rio.p- 



Proof. Part of the proof follows the same lines as for Theorem |3.14| We highlight here only the new elements 
needed for this proof. 



min n(z] 



V0l//(r2o,r) 



rain p{z) 



First, note that, for arbitrary r > and zq S 2?, 

(3.9) VOlMi^zo.r) ^ fJ'{z)dvo\Hiz) >VO\Hi^zo,r 

This motivates the definition of the sets OA,r, 

(3.10) Oa,- = \zeV \ nun p{z') > ^ 

[ z'en^^r volfl'(0, iZo.rJ 

A > is still arbitrary at this point; its value will be set below. 

Now pick r < R, and set e = (i? — r)/2. Note that if z G OA,r, then vo\m{^z,r) > vol7vi(f^2,r) > A. 

Since p is bounded below by a strictly positive constant on each rio,p'j we can pick, for arbitrarily large p, 
A > such that f2o,p C Oa.t: for this it suffices that A exceed a threshold depending on p and r. (Since 
p{z) — > as z approaches the boundary of V in Euclidean norm, we expect this threshold to tend towards 
as p — >■ oo.) We assume that fio.p C OA,r in what follows. 



Similar to the proof of Theorem 3.14 we invoke Lemma 3.13 to infer the existence of zq^wq such that 
zo € ^lo,e/2 and djj j^{zo,wo) — 0. We denote 

p* = sup{r' ; „(zo, Wo) = 0}; 

as before, there exists a Mobius transformation m such that h'{m{z)) = p{z) for all z in Clza,p''- To complete 
our proof it therefore suffices to show that p* > p + e/2, since Hq p C ^zo.p+s/2 ■ 



Suppose the opposite is true, i.e. p* < p + e/2. By the same arguments as in the proof of Theorem 3.14 
there exists, for each zi S Sfi^ 



n* — r — ^1 a 



point Z2 G ^zi.e/2 such that djj ^{z2,W2) — for some W2. Since 
the hyperbolic distance between Z2 and is bounded above by e/2 + p* — r — e + e/2 < p — r + e/2 < p, 
Z2 € f^o.p C OA,r, SO that V0I7V1 {^z2.r) > ^- It then follows from the conditions on Ai and Af that h'{m(z)) — 
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n{z) for all z in ^lz„.p* U ^Z2,R ^ ^zo.p* ^ ^zi,r+3s/2- Repeating the argument for all zi G dftzg^p*-r-e shows 
that h'{m{z)) — ^(z) can be extended to all z £ f^zo,p'+e/2, leading to a contradiction that completes the 
proof. □ 



4. Discretization and implementation 



To transform the theoretical framework constructed in the preceding sections into an algorithm, we need to 



discretize the relevant continuous objects. Our general plan is to recast the transportation eq. (3.7) as a 
linear programming problem between discrete measures. This requires two approximation steps: 

1) approximating the surface's Uniformization, and 

2) discretizing the resulting continuous measures and finding the optimal transport between the discrete 
measures. 

To show how we do this, we first review a few basic notions such as the representation of (approximations 
to) surfaces by faceted, piecewise flat approximations, called meshes^ and discrete conformal mappings; the 
conventions we describe here are the same as adopted in |18j . 



4.1. Meshes, mid-edge meshes, and discrete conformal mapping. Triangular (piecewise-linear) meshes 
are a popular choice for the definition of discrete versions of smooth surfaces. We shall denote a triangular 
mesh by the triple M = (V, E, F), where V = {wi}™ ^ C M.^ is the set of vertices, E = {cij} the set of edges, 
and F = {fi.j^k} the set of faces (oriented i ^ j ^ k). When dealing with a second surface, we shall denote 
its mesh by N . We assume our mesh is homeoniorphic to a disk. 

Next, we introduce "conformal mappings" of a mesh to the unit disk. Natural candidates for discrete 
conformal mappings are not immediately obvious. Since we are dealing with piecewise linear surfaces, it 
might seem natural to select a continuous linear maps that is piecewise affine, such that its restriction to 
each triangle is a similarity transformation. A priori, a similarity map from a triangular face to the disk has 
4 degrees of freedom; requiring that the image of each edge remain a shared part of the boundary of the 
images of the faces abutting the edge, and that the map be continuous when crossing this boundary, imposes 
4 constraints for each edge. This quick back of the envelope calculation thus allows 4|F| degrees of freedom 
for such a construction, with 4|i?| constraints. Since 3|F|/2 \E\ this problem is over constrained, and 
a construction along these lines is not possible. A different approach uses the notion of discrete harmonic 
and discrete conjugate harmonic functions due to Pinkall and Polthier ^23,i25j to define a discrete conformal 
mapping on the mid-edge mesh (to be defined shortly). This relaxes the problem to define a map via a 
similarity on each triangle that is continuous through only one point in each edge, namely the mid point. 
This procedure was employed in |18| : we will summarize it here; for additional implementation details we 
refer the interested reader (or programmer) to that paper, which includes a pseudo-code. 

The mid-edge mesh M = (V, E, F) of a given mesh M — {V,E,F) is defined as follows. For the vertices 
Vr G V, we pick the mid-points of the edges of the mesh M; we call these the mid-edge points of M. There is 
thus a e V corresponding to each edge e^ j G E. If and are the mid-points of edges in E that share a 
vertex in M, then there is an edge es,r G E that connects them. It follows that for each face fi.j^k £ F we can 
define a corresponding face ir,s,t € F, the vertices of which are the mid-edge points of (the edges of) fij^k', this 
face has the same orientation as fi,j,k- Note that the mid-edge mesh is not a manifold mesh, as illustrated 
by the mid-edge mesh in Figure [2j shown together with its "parent" mesh: in M each edge "belongs" to only 
one face F, as opposed to a manifold mesh, in which most edges (the edges on the boundary are exceptions) 
function as a hinge between two faces. This "lace" structure makes a mid-edge mesh more flexible: it turns 
out that it is possible to define a piecewise linear map that makes each face in F undergo a pure scaling (i.e. 
all its edges are shrunk or extended by the same factor) and that simultaneously flattens the whole mid-edge 
mesh. By extending this back to the original mesh, we thus obtain a map from each triangular face to a 
similar triangle in the plane; these individual similarities can be "knitted together" through the mid-edge 
points, which continue to coincide (unlike most of the vertices of the original triangles). 
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Mid-edge mesh zoom-in 
Figure 2. A mammahan tooth surface mesh, with the corresponding mid-edge mesh. 



To determine the fiattening map, we use the framework of discrete harmonic and conjugate harmonic func- 
tions, first defined and studied by PinkaU and Polthier [23, 25j in the context of discrete minimal surfaces. 
This framework was first adapted to the present context in |18j : this adaptation is explained in some detail 
in Appendix B. The flattening map is well-defined at the mid-edges V^. As shown in [T5] (see also Appendix 
B) the boundary of the mesh gets mapped onto a region with a straight horizontal slit (see Figure |3j where 
the boundary points are marked in red). We can assume, without loss of generality, that this slit coincides 
with the interval [—2,2] C C, since it would suffice to shift and scale the whole figure to make this happen. 
The holomorphic map z ^ w + ^ maps the unit disk conformally to C \ [—2, 2], with the boundary of the 
disk mapped to the slit at [—2, 2]; when the inverse of this map is applied to our flattened mid-edge mesh, 
its image will thus be a mid-edge mesh in the unit disk, with the boundary of the disk corresponding to the 
boundary of our (disk-like) surface. (See Figure [sj) We shall denote by <I> : V — ^ C the concatenation of 
these different conformal and discrete-conformal maps, from the original mid-edge mesh to the corresponding 
mid-edge mesh in the unit disk. 

Next, we define the Euclidean discrete conformal factors, defined as the density, w.r.t. the Euclidean metric, 
of the mid-edge triangles (faces), i.e. 



Note that according to this definition, we have 



vo1r3 (fr,s,t) 

vol($(f,.,.,t)) ■ 



E J , VOl]R3(fr,s.t) 

/i, dvolE ~ 



voIe («'(fr,5,t)) 



voIb (^(fr.s,*)) — volB_3(fr,s,t), 



where vol^ denotes the standard Euclidean volume element dx^ A dx^ in C, and volH3(f) stands for the area 
of f as induced by the standard Euclidean volume element in M^. The discrete Euclidean conformal factor 
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Mid-edge uniformization Uniformization Zoom-in 




After mapping to the disk Interpolated conformal factor 



Figure 3. The discrete conformal transform to the unit disk for the surface of Figure [2| 
and the interpolation of the corresponding discrete conformal factors (plotted with the JET 
color map in Matlab). The red points in the top row's images show the boundary points of 
the disk. 

at a mid-edge vertex is then defined as the average of the conformal factors for the two faces ir,s,t and 
ir,s',t' that touch in V^, i.e. 

Figure [3] illustrates the values of the Euclidean conformal factor for the mammalian tooth surface of earlier 
figures. The discrete hyperbolic conformal factors are defined according to the following equation, consistent 
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with the convention adopted in section [2] 

(4.1) ^il = (i-|$(v,,)P)'. 

As before, we shall often drop the superscript: unless otherwise stated, /i — /i^, and v = . 

The (approximately) conformal mapping of the original mesh to the disk is completed by constructing a 
smooth intcrpolant : I? — ^ K, which interpolates the discrete conformal factor so far defined only at the 
vertices in $(V); is constructed in the same way. In practice we use Thin-Plate Splines, i.e. functions of 
the type 

^^i{z) = Pi{z) +^b, ij{\z - z,\) , 

i 

where ''^{r) = log(r^), pi{z) is a linear polynomial in x^, a;^, and 6^ G C; pi and the bi are determined by the 
data that need to be interpolated. Similarly Tjj{w) = qi{w) + J^j '^j ''Pi\w ~ Wjl) for some constants Cj G C 
and a linear polynomial qi{w) in y^,y^. We use as interpolation centers two point sets Z — {ziY^^^ , and 
W = {wj}^^^ defined in the next subsection for the discretization of measures. See Figure js] (bottom-right) 
the interpolated conformal factor based on the black point set. 

We also note that for practical purposes it is sometimes advantageous to use Smoothing Thin-Plate Splines: 
r,(.) = argmin l.v„ - ,iH^>r))f + (1 - A) {^)\ A ..^} . 

when using these, we picked the value 0.99 for the smoothing factor A. 



4.2. Discretizing continuous measures and their transport. In this subsection we indicate how to 
construct discrete approximations T|?s(,j. (j(CiC) for the distance d^(X,y) = T^(^X) between two surfaces 
X and 3^, each characterized by a corresponding smooth density on the unit disk T) (f for X, ( for y). (In 
practice, we will use the smooth functions F^ and F^ for ^ and ) We shall use discrete optimal transport 
to construct our approximation T^^^^j. ^(^, C)i based on sampling sets for the surfaces, with convergence to 
the continuous distance as the sampling is refined. 

To quantify how fine a sampling set Z is, we use the notion oi fill distance ip{Z): 

ip{Z) := sup {r > I z e 7W : Bg{z, r) n Z,, = 0} , 

where Bg{z,r) is the geodesic open ball of radius r centered at z. That is, (p{Z) is the radius of the largest 
geodesic ball that can be fitted on the surface X without including any point of Z. The smaller (p{Z), the 
finer the sampling set. 

Given the smooth density ^ (on D), we discretize it by first distributing n points Z = {zi}2^i on X with 
ip{Z) — h > 0. For i — 1,. .. ,n, we define the sets to be the Voronoi cells corresponding to Zi e Z; 
this gives a partition of the surface X into disjoint convex sets, X — Uf^j^S^. We next define the discrete 
measure as a superposition of point measures localized in the points of Z, with weights given by the areas 
of Si, i.e. £,z = I]r=i C»<5zi, with := ^(S^) = dvol^- Similarly we denote by = {wj}^^^, T^, and 
(j := C(Tj) the corresponding quantities for surface 3^. We shall always assume that the surfaces X and y 
have the same area, which, for convenience, we can take to be 1. It then follows that the discrete measures 
£,z and (w have equal total mass (regardless of whether n — p oi not). The approximation algorithm will 
compute optimal transport for the discrete measures £_z and (^vk; the corresponding discrete approximation 
to the distance between ^ and ( is then given by T^{£,z, Cw)- 

Convergence of the discrete approximations T^i^zXw) to T^i^X) = ti^(A',3^) as ^{Z), ip{W) — t- then 
follows from the results proved in [17]. Corollary 3.3 in fTT requires that the distance function d^^{-, •) used 
to define T^{^, C) be uniformly continuous in its two arguments. We can establish this in our present case by 



invoking the continuity properties of d^\- proved in Theorem 3.7 extended by the following lemma, proved 
in Appendix A. 
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Lemma 4.1. Let {(z/t, Wk)}f^yi d V xV be a sequence that converges, in the Euclidean norm, to some point 
in {z' , w') £ V X 'D\'D X V, that is \zk — z'| + \wk — ui'j — >■ 0, as k ^ oo. Then, lim/^^ao d^^{zi.,Wk) exists 
and depends only on the limit point {z',w'). 

We shall denote this continuous extension oi dj}j^{-,-) to V x T) by the same symbol d^^. 

Since T) x D is compact, (this extension of) d^'^{-, •) is uniformly continuous: for all e > 0, there exists a 
6 = (5(e) such that, for all z, z' € X, w, w' e y, 

dxiz, z') < S{e) ,dy{w,w') < 5{e) ^ \dfj-{z,w) - dfj^{z' ,w')\ < e, 

where dx{-,-) is the geodesic distance on X, and dy{-,-) is the geodesic distance on y. 

The results in [T7] then imply that S^z ^ £, ^"^^ the weak sense, as (p{Z) 0, i.e. that for all bounded 
continuous functions / : I? ^> M, the convergence / d£,z — > /p / rf^ holds [221 . Similarly and Cw C in 
the weak sense as ip{W) 0. Furthermore, [17] also proves that for max{Lp{Z), (p{W)) < 

\T^i^zXw)-Tf{iX)\<s. 

More generally, it is shown that 

(4.2) K{^z,Cw)-T^{^,C)\<^dl^ (max(^(Z),^(W^))), 

where w^r is the modulus of continuity of d^^, that is 



5(£) 



sup 

d;f;(z,z') + dy(w,w')<t 



\dl^{z,w) 



dl^{z',w')\ 



We shall see below that it will be particularly useful to choose the centers in Z = {^iliLi: ^ — such 
that the corresponding Voronoi cells are (approximately) of equal area, i.e. n = N ~ p and S^i — ^(S^) « 
Cj = CC^i) ~ where we have used that the total area of each surface is normalized to 1. An effective way 
to calculate such sample sets Z and W is to start from an initial random seed (which will not be included 
in the set), and take the geodesic point furthest from the seed as the initial point of the sample set. One 
then keeps repeating this procedure, selecting at each iteration the point that lies at the furthest geodesic 
distance from the set of points already selected. This algorithm is known as the Farthest Point Algorithm 
(FPS) [7]. An example of the output of this algorithm, using geodesic distances on a disk-type surface, is 
shown in Figure [4] Further discussion of practical aspects of Voronoi sampling of a surface can be found in 



Figure 4. Samphng of the sur- 
face of Figure |2] obtained by the 
Farthest Point Algorithm. 



4.3. Approximating the local distance function d^^. We are now ready to construct our discrete 



version of the optimal volume transportation for surfaces (3.7). The previous subsection describes how to 



derive the discrete measures yiz^vw from the approximate conformal densities r^,r,y and the sampling 
sets Z and W . For simplicity, we will, with some abuse of notation, identify the approximations rp,r^ 
with /z, V. The approximation error made here is typically much smaller than the errors made in further 
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steps (see below) and we shall neglect it. The final component is approximating ^{zi,Wj) for all pairs 
[zijWj) Cz Z X W. Applying (3.4) to the points Zi, wj we have: 



(4.3) 



fi{z) — h'{m{z)) 



To obtain d^\,{zi,Wj) we will thus need to compute integrals over hyperbolic disks of radius R, which is 
done via a separate approximation procedure, set up once and for all in a preprocessing step at the start of 
the algorithm. 

By using a Mobius transformation fh such that ?7i(0) = Zq, and the identity 



^{m{u)) — ^{m o m(u)) dvolniu) = 



zo.R 



^{z)) — v{m{z)) dvolniz) 



we can reduce the integrals over the hyperbolic disks ^Zi,R to integrals over a hyperbolic disk centered around 
zero. 

In order to (approximately) compute integrals over Qq — ^o,/? — {z\ \z\ < vr}, we first pick a positive integer 
K and distribute centers pk-, k = 1,...,K in i7o- We then decompose Qq into Voronoi cells corresponding 
to the Pk, obtaining Hq = U^^^A^; see Figure [s] (note that these Voronoi cells are completely independent 
of those used in 4.2 ) 

To approximate the integral of a continuous function / over Qq we then use 

/ fiz)dvo\Hiz) / dvolif(z) f{pk)^^akf{pk) 

where au = /^^ dm\H(z). 

We thus have the following approximation: 



mm 



mm 



fi(z) ~ iy{m{z)) dvolniz) 
li{mi{z)) - vimlrhilz))) dmlniz) 



(4.4) 



min Qffc fi{mi{pk)) - v{m{mi{pk))) 
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where the Mobius transformations rrii, mapping to Zi, are selected as soon as the Zi themselves have been 
picked, and remain the same throughout the remainder of the algorithm. 

It can be shown that picking a set of centers {pk} with fill-distance h > leads to an 0{h) approximation; 
in Appendix ] A| we prove: 

Theorem 4.2. For continuously dijferentiable /i, ly, 



d^^{z^,Wj)- min | fi{mi{pk)) - i^{m{m^{pk))) \ 

m(zi)=Wj ^ 



<C^{{Pk}) 



where the constant C depends only on R. 



Let us denote this approximation by 



min Ok \ fJ,{mi{pk)) - i^{m{m^{pk))) 

m{zi)—Wj — 



Since the above theorem guarantees that the approximation error d^ ^{zi,Wj) — d^ ^{zi,Wj) 



can be uni- 



formly bounded independently of Zi , Wj , it can be shown that 

Tf{fiz,i^w)-T^il^z,^w)\<C^i{pk}), 
where again C is dependent only upon /i, v, R. Combining this with eq.( |4.2[ ) we get that 
(4.5) Tii^^,,y)~T^{^iz,l^w) < (max (^(Z), ^(T4^))) + ({pj) . 



In practice, for calculating d^^, the minimization over Mu zi^wj, the set of all Mobius transformations that 
map Zi to w/c, is discretized as well: instead of minimizing over all mz^^wj.a (see subsection 3.1), we minimize 
over only the Mobius transformations {razi^wj,2-Ki/L) ^ Taking this into account as well, we have 

thus 



(4.6) 



mm Ok 
'=i,...L ^ 

k 



lJ.{mi{pk)) - v{rnz 



Ti/Li-niiipk))) 



the error made in approximation (4.6 1 is therefore proportional to L -'^ -I- Cip {{pk})- 



To summarize, our approximation T~{iJ,z,i^w) to the uniformly continuous T^{^,v) is based on two ap- 
proximations: on the one hand, we compute the transportation cost between the discrete measures ^z, i^w, 
approximating /x, i/; on the other hand, this transportation cost involves a local distance which is it- 
self an approximation. The transportation between the discrete measures will be computed by solving 
a linear programming optimization, as explained in detail in the next subsection. The final approxima- 



tion error (4.5) depends on two factors: 1) the fill distances (p{Z),ip{W) of the sample sets Z,W, and 
2) the approximation of the local distance function d^ ^{zi^Wj) between the sample points. Combining 



the discretization of the Mobius search with (4.6), the total approximation error is thus proportional to 

Recall that we are in fact using F^, in the role of of /i, v (see above), which entails an additional approxi- 
mation error. This error relates to the accuracy with which discrete meshes approximate smooth manifolds, 
as well as the method used to approximate uniformization. We come back to this question in Appendix B. 
As far as we are aware, a full convergence result for (any) discrete uniformization is still unknown; in any 
case, we expect this error to be negligible (and approximately of the order of the largest edge in the full 
mesh) compared to the others. 
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4.4. Optimization via linear programming. The discrete formulation of eq. (3.7) is commonly formu- 
lated as follows: 



(4.7) 



dijTTij min 

hi 



(4.8) 



TTin > 



where — and i^j — iy{Tj), and dij — d^jj[zi,Wj). 

In practice, surfaces are often only partially isometric (with a large overlapping part), or the sampled points 
may not have a good one-to-one and onto correspondence (i.e. there are points both in Z and in W that 
do not correspond well to any point in the other set). In these cases it is desirable to allow the algorithm 
to consider transportation plans tt with marginals smaller or equal to /i and v. Intuitively this means that 
we allow that only some fraction of the mass is transported and that the remainder can be "thrown away" . 
This leads to the following formulation: 



(4.9) 



(4.10) 



TT y > 



where < Q < 1 is a parameter set by the user that indicates how much mass must be transported, in total. 

The corresponding transportation distance is defined by 

(4.11) Tdiv, J^) = X! ^y '^y ' 



where iTij are the entries in the matrix tt for the optimal (discrete) transportation plan. 
Since these equations and constraints are all linear, we have the following theorem: 



Theorem 4.3. The equations (4.7)-(4.8) and (4.9)- (4.10) admit a global minimizer that can be computed 



in polynomial time, using standard linear-programming techniques. 



When correspondences between surfaces are sought, i.e. when one imagines one surface as being transformed 
into the other, one is interested in restricting tt to the class of permutation matrices instead of allowing all 
bistochastic matrices. (This means that each entry ir^j is either or 1.) In this case the number of centers 



Zi must equal that of wj, i.e. 



N — p, and it is best to pick the centers so that /i^ 



N 



Vj, for all 



i, j. It turns out that this is sufficient to guarantee (without restricting the choice of tt in any way) that the 
minimizing tt is a permutation: 



Theorem 4.4. If n — N ^ p and 



N 



then 



(1) There exists a global minimizer o/(4.7) that is a permutation matrix. 



(2) If furthermore Q = where M < N is an integer, then there exists a global minimizer of iA.9h tt 



such that TTij G {0, 1} for each i, j. 



Remark 4.5. In the second case, where TTij G {0, 1} for each i, j and ~ ^ '■^^^ ^^i^l t)e viewed 

as a permutation of M objects, "filled up with zeros" . That is, if the zero rows and columns of tt (which 
must exist, by the pigeon hole principle) are removed, then the remaining M x M matrix is a permutation. 
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Proof. We first note that in both cases, we can simply renormalize each fii and i^j by N, leading to the 
rescaled systems 



(4.12) 




E 

E 



3 -^3 

n, 



> 



< 1 

< 1 

= M 



To prove the first part, we note that the left system in (4.12) defines a convex polytope in the vector space 
of matrices that is exactly the Birkhoff polytope of bistochastic matrices. By the Birkhoff-Von Neumann 
Theorem |15) every bistochastic matrix is a convex combination of the permutation matrices, i.e. each tt 



satisfying the left system in (4.12) must be of the form Efc'^^''''^' where the are the A^! permutation 



matrices for N objects, and Ck — 1, with Ck > 0. The minimizing tt in this polytope for the linear 
functional (4.7) must thus be of this form as well. It follows that at least one t'' must also minimize (4.7 1, 



since otherwise we would obtain the contradiction 



(4.13) 



> mm 

fc 



> 



The second part can be proved along similar steps: the right system in (4.12) defines a convex polytope 



in the vector space of matrices; it follows that every matrix that satisfies the system of constraints is a 
convex combination of the extremal points of this polytope. It suffices to prove that these extreme points are 
exactly those matrices that satisfy the constraints and have entries that are either or 1 (this is the analog 
of the Birkhoff-von Neumann theorem for this case; we prove this generalization in a lemma in Appendix 
C); the same argument as above then shows that there must be at least one extremal point where the linear 
functional (4.7) attains its minimum. □ 



This means that, when we seek correspondences between two surfaces, there is no need to impose the (very 
nonlinear) constraint on n that it be a permutation matrix; one can simply use a linear program to solve 
either , with Theorem 4.4 guaranteeing that the minimizer for the "relaxed" problem (4.7)-(4.8) or (4.9)- 
(4.10) is of the desired type if n = N — p and fii = = vj. 



4.5. Consistency. In our schemes to compute the surface transportation distance, for example by solving 



(4.9), we have so far not included any constraints on the regularity of the resulting optimal transportation 
plan TT* . When computing the distance between a surface and a reasonable deformation of the same surface, 
one does indeed find, in practice, that the minimizing tt* is fairly smooth, because neighboring points have 
similar neighborhoods. There is no guarantee, however, that this has to happen. Moreover, we will be 
interested in comparing surfaces that are far from (almost) isometric, given by noisy datasets. Under such 
circumstances, the minimizing tt* may well "jump around" . In this subsection we propose a regularization 
procedure to avoid such behavior. 

Computing how two surfaces best correspond makes use of the values of the "distances in similarity" 
^{zi,Wj) between pairs of points that "start" on one surface and "end" on the other; computing these 
values relies on finding a minimizing Mobius transformation for the functional (3.4). We can keep track of 
these minimizing Mobius transformations niij for the pairs of points {zi,Wj) proposed for optimal correspon- 
dence by the optimal transport algorithm described above. Correspondence pairs (i, j) that truly participate 
in some close-to-isometry map will typically have Mobius transformations rriij that are very similar. This 
suggests a method of filtering out possibly mismatched pairs, by retaining only the set of correspondences 
that cluster together within the Mobius group. 

There exist many ways to find clusters. In our applications, we gauge how far each Mobius transformation 



is from the others by computing a type of ii variance: 



(4.14) 



Evi^,j) 



E 

(k,e) 



mkel 
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where the norm is the Frobenius norm (also caUed the Hilbert-Schmidt norm) of the 2x2 complex matrices 
representing the Mobius transformations, after normalizing them to have determinant one. We then use 
Ey{i,j) as a consistency measure of the corresponding pair (i, j). 

5. Examples and comments 






(a) 



Good pair (a) 





(b) 



Erroneous pair (b) 



Figure 6. Calculation of the local distance dj}^{-,-) between pairs of points on two different 
surfaces (each row shows a different pair of points; the two surfaces are the same in the top 
and bottom rows) . The first row shows a "good" pair of points together with the alignment 
of the conformal densities ^^wfu based on the best Mobius transformation m minimizing 
11^ — 771*1^11 dvolTvi. The plot of this latter integral as a function of m (parameterized by 



a G [0, 27r), see (2.3)) is shown in the right-most column. The second row shows a "bad" 
correspondence which indeed leads to a higher local distance ^. 

In this section we present a few experimental results using our new surface comparison operator. These 
concern an application to biology; in a case study of the use of our approach to the characterization of 
mammals by the surfaces of their molars, we compare high resolution scans of the masticating surfaces of 
molars of several lemurs, which are small primates living in Madagascar. Traditionally, biologists specializ- 
ing in this area carefully determine landmarks on the tooth surfaces, and measure characteristic distances 
and angles involving these landmarks. A first stage of comparing different tooth surfaces is to identify 
correspondences between landmarks. Figure [g] illustrates how d^j^{z,w) can be used to find corresponding 
pairs of points on two surfaces by showing both a "good" and a "bad" corresponding pair. The left two 
columns of the figure show the pair of points in each case; the two middle columns show the best fit af- 
ter applying the minimizing Mobius on the corresponding disk representations; the rightmost column plots 
I fi{z) — irnl^^^^^„v){z) I dvolniz), the value of the "error" , as a function of parameter u, parametrizing 



the Mobius transformations that map a give point zq to another given point wq (see Lemma 3.5 1. The "best 



corresponding point wq for a given zq is the one that produces the lowest minimal value for the error, i.e. 
the lowest dj},^(zQ,wo). 

Figure [t] show the top 120 most consistent corresponding pairs (in groups of 20) for two molars belonging 
to lemurs of different species. Corresponding pairs are indicated by highlighted points of the same color. 
These correspondences have surprised the biologists from whom we obtained the data sets; their experimental 
measuring work, which incorporates finely balanced judgment calls, had defied earlier automatizing attempts. 
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Once the differences and similarities between molars from different animals have been quantified, they can 
be used (as part of an approach) to classify the different individuals. Figure |8] illustrates a preliminary 
result from [T3] that illustrates the possibility of such classifications based on the distance operator between 
surfaces introduced in this paper. The figure illustrates the pairwise distance matrix for eight molars, coming 
from individuals in four different species (indicated by color) . The clustering was based on only the distances 
between the molar surfaces; it clearly agrees with the clustering by species, as communicated to us by the 
biologists from whom we obtained the data sets. 

One final comment regarding the computational complexity of our method. There are two main parts: the 
preparation of the distance matrix dij and the linear programming optimization. For the linear programming 
part we used a Matlab interior point implementation with unknowns, where N is the number of points 
spread on the surfaces. In our experiments, the optimization typically terminated after 15 — 20 iterations 
for N = 150 — 200 points, which took about 2-3 seconds. The computation of the similarity distance dij 
took longer, and was the bottleneck in our experiments. If we spread N points on each surface, and use 
them all (which was usually not necessary) to interpolate the conformal factors r^,r^, if we use P points 
in the integration rule, and take L po ints in the Mobius discretization (see Section |4] for details) then each 



approximation of d^,^{zi,Wj) by (4.6) requires 0{L ■ P ■ N) calculations, as each evaluation of T^,Ti, takes 
0{N) and we need L ■ P oi those. Since we have 0{N'^) distances to compute, the computation complexity 
for calculating the similarity distance matrix dij is 0{L ■ P ■ N^). In practice this step was the most time 
consuming and took around two hours for TV = 300. However, we have not used any code optimization and 
we believe these times can be reduced significantly. 
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Appendix A. 



This Appendix contains some technical proofs of Lemmas and Theorems stated in sectio n [3 



with proving the list of properties of the distance function d^^{z,w) given in Theorem 
Theorem 13.3 



3.3 



and|4] We start 
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The distance function d^'^(z,w) satisfies the following properties 

(1) d^if_t,m*,^i''^'i^i^o),'m'2^{'Wo)) = d^j^{zo,wo) Invariance under (well-defined) 

Mobius changes of coordinates 

(2) d^,„{zo,Wo) ^ d^'^^{wo,zo) Symmetry 

(3) d^^,{zn,wo)>0 Non-negativity 

(4) dj}j^{zo,wo) = ill ^-^id i^wo,R in [T^^v] a-i'e isometric 

(5) c^m•l..l/('7^~^(^o),2;o) = Refiexivity 

(6) <,M3( 

^ii 23) < dj}_^^^^{zi, Z2) + dj}^^^^{z2, 23) Triangle inequality 



Proof. For (1), denote "'^(^o) = ^1, and m2 ^{wq) — Wi. Then 

'^mTM,™>(^i''^i) = ^i'^f / |mj^(z) - m*TO;zy(z)|dvol/f(z) 

[if / |^(toi(z)) - j/(m2(m(z)))|dvol/f-(z). 



: inf 

m{z-i 



Next set m = m2 o mo nii ^. Note that m(zo) = wo- Plugging TO2(m(z)) = m(mi(z)) into the integral and 
carrying out the change of variables mi(z) ~ z' , we obtain 

inf / l^(z') - i^(TO(z')) Mvol^(z') = ^ inf / \[i{z') - v{m{z'))\d^o\H{z'). 



For (2), we use Lemma 3.2 and equations (2.5), (2.6) to write 

d^^{zo,Wo)^ inf / |/i(z) - mV(z)|dvol/f(z) 

= inf / \{Tn~^)*n{w)-iy{w)\dvolH{w)^d^{wo,zo). 



(3) and (4) are immediate from the definition of d^^. 



(5) follows from the observation that the minimizing m (in the definition (3.4) of dj} ^) is nii itself, for which 
the integrand, and thus the whole integral vanishes identically. 

For (6), let mi be a Mobius transformation such that mi(zi) = Z2, and m2 such that 7712(22) = -^s- Setting 
m — TO2 ° "^1 1 we have 



(A.l) 



d^,^f,^izi,Z3) < / I /^l(z) - TO>3(2) I rfvolH(z) 



i.i? 



< / I /xi(z) - mi^2(z) I rfvoli/(z) + / \mlfi2{z) ~m*fj,3{z)\dmlH{z) . 
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The second term in ( A.l I can be rewritten as (using Lemma 3.2 the change of coordinates r7ii(zi) — zi and 
the observation m* = 77117712) 



7nlii2{z) — 7n* ii:i{z) \ dvolniz) — / | TOi*m2/i2(u;) — mi*TO^m2/i3(w) | (ivol//(w) 

,-R ^^Z2,R 

I ^i2iw) - 7n*2^l3{w) I dvolH{w). 



Z2.R 

We have thus 

f^Ml,M3(^l'^3) - / \fJ.i{z) -7nlfJ.2{z)\dvolHiz) + \ fl2iw) ~ 7n*2^3{w) \ dYolniw) , 

and this for any mi, 7712 € Md such that r7ii(zi) = Z2 and 7712(2:2) = 23. Minimizing over ttii and 7712 then 
leads to the desired result. 

□ 

Next we prove the continuity properties of the function $(zo, w^o, o") = I /^(-z)— ^^(™zo,too,(T(2:)) | dvoljy (z), 

stated in Lemma 3.6 which were used to prove continuity of d^ ^ itself (in Theorem 3.7). 

Lemma 13.61 

• For each fixed {zq,wq) the function $(zoj Wo, ') is continuous on Si. 

• For each fixed a E Si, $(•, •, ct) is continuous on D x D. Moreover, the family ( $(•, •, cr) J is equicon- 
tinuous. 

Proof. We start with the continuity in a. We have 

\^{zo,wo,a-) - ^{zo,wo,a') \ < / 11^(777^0. ^(z)) - ;/(7772o.^o,cr' (2)) M™I/f-(z) . 

Jnizo,R) 



Because v is continuous on T), its restriction to the compact set fl{wo, R) (the closure of fl{wo, R)) is bounded. 
Since the hyperbolic volume of ^{zq, R) is finite, the integrand is dominated, uniformly in cr', by an integrable 
function. Since mzg.wo.a-i^) is obviously continuous in cr, we can use the dominated convergence theorem to 
conclude. 

Since S^ is compact, this continuity implies that the infimum in the definition of djj\, can be replaced by a 
minimum: 

d^^{za,wa)^ min / \ ^i{z) - v{m{z))\dvo\H(z) . 

Next we prove continuity in zq and wq (with estimates that are uniform in cr). 

Consider two pairs of points, {zq,wq) and {z'q,w'q) eV x V. Then 
\<i>{zo,WQ,cr) ~'i>{z'o,w'o,a)\ 



I ^(z) - 7/(?77^„^^„,cr(z)) I dvol//(z) - / \ ^{u) - iy{m,,'^^^'^^{u) \ dvolH (u) 

n(zo..R) Jn(z'g,R) 

< I |A*(^)-'^("lzo.«'0,T(2))|rfvol//(z)- / |^(7772„^^J^iz) - I/(7772._„^_^ 7772g_2j^l(z)) |dvoljf(?7) 

0,-R) Jn{zo,R) 

< I (I - M(w^o,4,i(^)) I + I '^{mzo,wo,a{z)) - J^(77l3^.^J,.a(m^„.z^,,l(^))) I ) rfvola(z) . 
JQ{zo,R) 

On the other hand, note that for any 7 > 0, /i and v are continuous on the closures of fl{zQ,R + 7) and 
i}{wo, R + j), respectively; since these closed hyperbolic disks are compact, /i and v are bounded on these 
sets. Pick now p > such that |zg — zo| < p, \w'q — wo\ < p imply that ri(zQ, R) C f^(zo, i? + 7) as well as 
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f2(u'Q, R) C i^{wQ, R + 7). It follows that, if \zq — zq\ < p and \w'q — wq\ < p, then | — p{mza,z'^,i(,z) \ and 
I ^ij^zoo,wo.a{z)) — i^{^z[^^w^<y{^zo,z'g,i{z))) I ^i'^ bounded uniformly for z £ fl{zQ,R). Since it is clear from 
the explicit expressions (3.6) that rnzg^z^.iiz) — ^ z and iTiz^^w'g,cri'm'zo,z'g.i{z)) — ^ mzg.wa.criz) as Zq — ^ zq and 
w'q — >■ Wo, we can thus invoke the dominated convergence theorem again to prove continuity of $(•, •,(t). 

To prove the equicontinuity, we first note that v is uniformly continuous on fl(wo,R) U fl{wQ,R), since v 
is continuous on the compact set ^l{wo, R + j), which contains il{wo,R) U J7(wQ,i?) for all w'q that satisfy 
\wq — wo\ < p. This means that, given any e > 0, we can find ^ > such that — i^{w')\ < e holds for 

all w, w' that satisfy w, w' € U,{wo^ R) U ^(w'q, R) and |w — w'j < 6. This implies the desired equicontinuity 
if we can show that Imzo^wo.criz) — 'mz'g,w'o,<TiiTT'zo,z' ,iiz))\ can be made smaller than (5, uniformly in cr G 5*1, 
by making |zq — zqI + |wq — WqI sufficiently small. 

We first estimate \mzo,wo,a{z) — mz„^w;-,^aiz)\- With the notations of (3.6), we have 

, (zq - WQa){l -z^w'^a) - {zq - WQa){l -z^WoW) 

a{zo,WQ,a) " a(zo,Wo,cr) = y- — — — _ , _ 

(1 - zoWoa)(l - zow'^a) 

^ (wq - WqMIzqP - 1) 
{1 -ZoWoW){l -Z^WqW) 

so that 

. , \ / / M I ^0 — W'n I ^ 

\a{zo,wo,a)-a{zo,Wo,a)\ < j- — p r— —r — — y < — p ^ —r — -y 

(1 - |zo| \wo\)[l - \zo\{\wo\ + Ol (1 - kol \wo\)[l - l^oKhol + 0] 

when |wo — w'q\ < It thus suffices to choose so that < ({1 — \zo\ |wo|)[l — koKl^^ol + 0] to ensure that 
\a{zQ,wo,a) ~ a(zo, Wo,cr)| < C- For the phase factor r in (3.6) we obtain 

(1 - zoWQa){l - zow^a) - (1 - zoWoa){l - ZQW^a) 



T{zo,wo,a) - T{zo,WQ,a) = a 



= a 



(1 - zoWoct)(1 - zoWqO-) 

(wq - w'q)z^- (w^- w[,)zog + \zq\^{w^w'q - w'qWq) 

(1 - z^woa){l - z^w'^^a) 
{m - w'q)z^ - zn{w^ - w'i^)a + \zo\'^[w^{w'g - Wq) + Wo{w^ - w'^]) 



(1 - zowoo-)(l - zoWqct) 
when Itdq — < ^, this implies 

I . ^ , ' |zo||u;n|[2+|zo|(2|^;o|+0] ^ 

(1 - |zo| |wo|)[l - \zq\{\wo\ + £,)] 
which can clearly be made smaller than any C > by choosing ^ sufficiently small. All this implies that (use 
^) 

l + \z\ (l + |z|)2 

\'mzo,wo,a{z) ~ rnzg^ni'g,a{z)\ < \T{zo,Wo,a) - t(zo, w[),cr)| ^ _ + \a{zo,Wo,a) - a(zo, Wq, cr)| ^j— 

^ . 2(l + |z|) 
(l-|z|)2' 

which will be smaller than S/2, uniformly in tr, if ^ < 6{1 — |zp)/8; this bound on in turn determines 
the bound to be imposed on the ^ used above. Hence Imzg^wo.aiz) — rnzg.w'g.aiz)\ < S/2 can be guaranteed, 
uniformly in cr, by choosing |wo — Wq| < ^ for sufficiently small ^. 
One can estimate likewise 

|™^0,<,t(^) - mz'^,W'gA^Z'g,ZO,liz))\ , 

and show that this too can be made smaller than S/2, uniformly in ct, by imposing sufficiently tight bounds 
on |zq — zqI and jwQ — wol. Combining all these estimates then leads to the desired equicontinuity, as indicated 
earlier. □ 



To prove Lemma |4.1[ we shall use the following lemma: 
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Lemma A.l. Consider Uk ~ e''^ + Sk, where |£/c| as fc ^ oo. Then there exists, for every e > 0, a 
-fsT £ IN such that for all k > K ; and all m G -A/^i.o.ufc , 

inf |m(w)| > 1 — £. 



The set Mjj^a.uk used in this lemma is given by Definition 3.4 



Proof. From Lemma 3.5 we can write m as 



fhiw) = e 



-10 



1 + UkC^^W ' 

for some G [0, 2tt). Substituting — e^"^ + in this equation we get 

m[w) — e —e^ ^ — -. — 

l + {e^^ +ek)e^<^w I + we^^<^-^) +eie' 

Writing the shorthand s for s = 1 have thus 



\0, 



rh{w) — e 



i4> 



s + sTkC^^w 



< 



< 



Eke-'^f' - Ske'^w 



\s + efce'^w 



< 



\ek\{l + \w\) 



\£k\\w\ 



Now for aU w e ^o,R, \w\ < rji ~ tanh ^{R). This imphes |s| > 1 — |w| > 1 — rji, and 1 + |w| < 1 + r_R, so 
that 



fh{w)-e'^ <\ek\ 



'l-rR-\ek\rR '^''^ I - rn {1 + \ek\)' 



Since \£k\ the lemma follows. 



□ 



We are now ready for 

Lemma |4.1 Let {{zk,Wk)}i,yi (I'D xV be a sequence that converges, in the Euclidean norm, to some point 



(z', w') G D X 'D\'D X D, that is \zk — z'| + \wk — w'| — > 0, as k ^ oo. Then, limfc_>.oo df,^{zki Wk) exists and 
depends only on the limit point [z' ,w'). 



Proof. Since (z', w') gV xV\V xV either z' gV\V oi w' gV\V. Let us assume that z' gV\V (the 
case w' € 2?\ I? is similar). Denote by an arbitrary Mobius transformation in M/j.o.iofc- By symmetry of 
the distance and using a change of variables we then obtain 



df(-{zk,Wk) = df^^{wk,Zk) 



mm 

m(wk)=Zk Jn 



dvolniw) 



mm 



C(mfc(w)) - ^(m(mfe(w))) dvolniw) 



Now, recall that ^(z) = £,^{z) = ^(z)(l — I^P)^, where £,{z) is a bounded function, sup^gp IC(^)I ^ C*?. From 



Lemma A.l we know that for every e > and for k > K sufhciently large, \m{mk{{w))\ > 1 — £ for all 
w S 1^0, ii J and all m such that m{wk) = Zfe- This means that for these k > K we have 



|2n2 



|^(m(mfe(w)))| = i{m{mk{w))) [l - \m{mk{w))\ ) 
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for all w e flo,R- Therefore, 



d?Azk,Wk)- I C{mk{wj) dvolniw) 



< 



< 



mm 



C(mfe(w)) - ^{m{mkiw))) dvolniw) - / C{mk{w)) dvolH{w) 



On 



mm / I ({rnkiw)) - £,{m{mk{w))) - C{mk{w)) \ dvolniw) 



< min 



^{■m{mk{w))) dvolniw) — > 0, as fc — )■ oo. 



Therefore d^i-{zk,Wk) converges, as A: ^ oo, if and only if J^^o r \C{^k{w))\ dvolH{w) converges, and to the 
same limit, for any mk G Mo.o.w^- We can take, for instance, mk{w) — which gives 



\C{m.k{w))\dwo\H{w) = 



W + Wk 



dm\fj(w). 



1 + WkW ^ 

For w G fio,i?,7 |1 + 'Wkw\ > 1 — rp. It follows that this expression has a limit as fc — )■ oo, and 

w + w' 



lim 



\C{mk{w))\ dvolniw) 



which clearly depends on w' , not on the sequence {wk}- 



1 + w'w 



dvolniw), 



□ 



Next, we prove Theorem |4.2[ We start with a simple lemma showing that all Mobius transformations 
restricted to Qo,r, R < oo, are Lipschitz with a universal constant, for which we provide an upper bound. 

Lemma A. 2. A Mobius transformation m G Mjj restricted to ^Iq^r, R < oo is Lipschitz continuous with 
Lipschitz constant Cm < (il^ilLl)'' ■ 



Proof. Denote m(z) = e j^r^- Then, for z,w <E ^o,R we have 



m{z) — m(w) 



< 



< 



-,16 



1 — za 1 — wa 

2 



< 



(z — a)(l — wa) — [w — a)(l — za) 



(1 — za){l — wa) 



iz-w)il-\an 
(1 — zd){l — wa) 



< \z — w\ 



il-rR\a\)^ 



□ 



Next we prove: 

Theorem |4.2| For continuously differentiable /x, v, 



min ak \ fJ.{mi{pk)) - v{m{mi{pk))) \ 

■{Zi)=Wj ^ 



<Cv{{Pk)) 



where the constant C depends only on ^, v, R. 

Proof. First, denote f{z) = ii(fhi{z)) — i/(TO(mi(z))) 
(A.2) 



Then, 



/ f{z)dvolH{z) - min Vafc/bfc) < Efc /o„ /(^) " /b^) '^^Iff (2) 
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where the modulus of continuity ll>^° (h) — sup|2_tu|</i.2 wecio ~ used. Note that 

(A.3) < uj"°~ ~ . 

V "/ / — fj.omi ' vomorrii 

Since /x, v have continuous derivatives on compact domain, they are Lipschitz continuous. Denote their 
Lipschitz constants by C^, C^, respectively. From Lemma A. 2 we see that, for z,w £ f^oi 

|2 



m^{z) - mi{w) 



1 



a-rR\a\y 



\z-w\< C, 



1 



which is independent of rrii . Similarly, 



j/(m(mi(z))) — v{m{'mi{'w))) 



< a 



m{mi(z)) — m(mi{w)) 



< a 



1 



I2: — w| 



which is independent of m,mi. Combining these with eq. ( A.2||A.3 ) we get 



f f{z)dvolH{z)^ min Vafc/(pfe) <{C^ + C,)jf 
which finishes the proof. 



2</'({Pfc}) : 



□ 



Appendix B. 



In this appendix we provide a short exposition on discrete and conjugate discrete harmonic functions on 
triangular meshes as presented in [6l 1231 124[ I25j , and we show how this theory can be used in our context to 
conformally flatten disk- type ( or even just simply connected) triangular meshes. 



We will use the same notations as in Section [4j Discrete harmonic functions are defined using a variational 
principle in the space of continuous piecewise linear functions defined over the mesh PLj^ as follows. 

Let us denote by (f>i{z), i = l,..,m, the scalar functions that satisfy 4>j{vi) = 6ij and are linear on each 
triangle /i S F. Then, the (linear) space of continuous piecewise-linear function on M can be written in 
this basis: 



PLm = \^Ui4i^{z) I (wi, ...,u„)'^ e 



Next, the following quadratic form is defined over PLm- 

(B.l) ED^riu) = J2[ (V^'V") 

where (•) = {-j-g^s denotes the inner-product induced by the ambient Euclidean space, and dvolgs is the 
induced volume element on /. This quadratic functional, the Dirichlet energy, can be written in coordinates 
of the basis defined earlier as follows: 



(B.2) 



J2 / (V0.,V</.,) 



M 



(V0i, V(t)j) dvolR3. 



The discrete harmonic functions are then defined as the functions u € PLm that are critical for Ejjir^u), 
subject to some constraints on the boundary of M. The linear equations for discrete harmonic function 
u G PLm are derived by partial derivatives of Eoir, (B.2) w.r.t. Ui,i = 1, ..,m: 



(B.3) ^^^^2Tu.. 



duk 



i=l 



^ 2 



= 2 



M 



where i?^ C M is the 1-ring neighborhood of vertex Vk- The last equality uses that (^^ is supported on i?fc. 

Now, let u — Ui(^i be a discrete harmonic function. Pinkall and Polthier observed that conjugating the 
piecewise-constant gradient field Vu (constant on each triangle f ^ F), i.e. rotating the gradient Vw in each 
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triangle / by 7r/2 in the positive ( = counterclockwise) sense (we assume M is orientable), results in a new 
vector field *du — Jdu with the special property that its integrals along (closed) paths that cross edges only 
at their mid-points are systematically zero (see for example [25] ). This means in particular that we can define 
a piecewise linear function *u such that its gradient satisfies d*u = *du and that is furthermore continuous 
through the mid-edges V G V. The space of piecewise-linear functions on meshes that are continuous through 
the mid-edges is well-known in the finite-element literature, where it is called ucPLm, the space of non- 



conforming finite elements [3j. The Dirichlet form (B.l ) is defined over the space of non-conforming elements 
ucPLm as well; the non-conforming discrete harmonic functions are defined to be the functions v G ucPLm 
that are critical for Ejjir and that satisfy some constraints on the mid-edges of the boundary of the mesh. 
Polthier [5S] shows that if m S PLm is a discrete harmonic function, then *u € ucPLm is also discrete 
harmonic, with the same Dirichlet energy, and vise-versa. Solving for the discrete harmonic function after 
fixing values at the boundaries amounts to solving a sparse linear system which is explicitly given in |25j . 

This theory can be used to define discrete conformal mappings, and used to flatten a mesh in a "discrete 
conformal" manner, as follows. The flattening is done by constructing a pair of conjugate piecewise linear 
functions (m, *m) where u G PLm, *u G ucPLm, and the flattening map $ : M — > C is given by 

(B.4) <i> = u + i * u. 

Since d * u = Jdu, $ is a similarity transformation on each triangle f £ F. Furthermore, $ is continuous 
through the mid-edges G V. This means that 3> is well-defined on the mid-edges V and maps them to the 
complex plane. 

The function u is defined by choosing an arbitrary triangle fout G F, excising it from the mesh, setting 
the values of u at two of /out's vertices Ui-^,Ui^ to and 1, respectively, and then solving for the discrete 
harmonic u that satisfies these constraints. See for example Figure [3] (top- left); the "missing mid-edge face" 
corresponding to the excised face fout would have connected the three mid-edge vertices that have a only 
one mid-edge face touching them. The conjugate function *u is constructed by a simple conjugation (and 
integration) process as described in [25] and [18]. 

A surprising property of the Discrete Uniformization <f> as it is defined above, which nicely imitates the 
continuous theory (see [22]) is that it takes the boundaries of M to horizontal slits, see Figure [3) top row 
(boundary vertices colored in red). This property allows us to easily construct a closed form analytic map 
(with "analytic" in its standard complex analytic sense) that will further bijectively map the entire complex 
plane C minus the slit to the open unit disk, completing our Uniformization procedure. 

This property is proved by arguments similar to those for Proposition 35 in j25] ; see also [18] . More precisely, 
we have 

Theorem B.l. Let $ : /W — > C be the flattening map from the mid-edge mesh M of a mesh M with bound- 
ary, using a discrete harmonic and conjugate harmonic pair as described above. Then, for each connected 
component of the boundary of A4, the mid-edge vertices of boundary edges are all mapped onto one line 
segment parallel to the real axis. 



Proof. Suppose u — '^iUi(j)i{-) is a discrete harmonic, piecewise linear and continuous function, defined at 
each vertex Vi G V , excluding the two vertices of the excised triangle for which values are prescribed; then 



we have, by (B.3 ) 



(B.5) 



(V(/)i, Vu)(ivolR3 ~ 0, 



Next, consider a boundary vertex Vj of the mesh Jv[. Denote by Vr,Vs the mid-edge vertices on the two 
boundary edges touching vertex Vj. We will show that *u(Vr) = *u(Vs); this will imply the theorem, since *u 
gives the imaginary coordinate for the images of the mid-edge vertices under the fiattening map (see (B.4)) . 
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Observe that on the triangle fi.j^k, 

(B.6) V</)j- = r. 

RecaUing that V * u = JVw, usmg |R6l, and j"^ = -J, we obtain 

= E '-V-^J- I/) 

= - / (Vw, V(/)j-) dvolR3 

0, 

where 7 is the piecewise hnear path starting at and passing through the mid-edge vertices of the 1-ring 
neighborhood of Vj ending at V^. The last equality is due to (B.3|. □ 

A natural question, when dealing with any type of finite-element approximation, concerns convergence as 
the mesh is refined: convergence in what sense, and at what rate? For discrete harmonic functions over 
meshes, this convergence is discussed in [121 El]. Note that these convergence results are in the weak sense; 
this motivated our defining the discrete conformal factors /Xf via integrated quantities (volumes) in Section 

12 

Finally, we note that the method presented here for Discrete Uniformization is just one option among 
several; other authors have suggested other techniques; for example [lOj . Typically, this part of the complete 
algorithm described in this paper could be viewed as a "black box" : the remainder of the algorithm would 
not change if one method of Discrete Uniformization is replaced by another. 



Appendix C. 



In this Appendix we prove a lemma used in the proof of Theorem 4.4 
Lemma The N x N matrices n satisfying 



(C.l) 



-Kij > 



constitute a convex polytope V of which the extremal points are exactly those tt that satisfy all these con- 
straints, and that have all entries equal to either or 1. 

Remark. Note that the matrices tt G V with all entries in {0, 1} have exactly M entries equal to 1, and all 
other entries equal to zero; if one removes from these matrices all rows and columns that consist of only 
zeros, what remains is a M x M permutation matrix. 



Proof. V can be considered as a subset of M , with all entries nonnegative, summing to M . The two 



inequalities in (C.l I imply that the entries of any n £ V are bounded by 1. These inequalities can also be 
rewritten as the constraint that every entry of AV — & e is non positive, where A is a 1 
and 6 is a vector in M^^. It follows that 7-" is a (bounded) convex polytope in . 



f2N 



matrix. 
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If TT e P C has entries equal to only or 1, then tt must be an extremal point of V by the following 
argument. If vr^ = 1, and tt is a nontrivial convex combination of tt^ and tt^ in 7^, then 

TT = A TT^ + (1 - A) TT^ with A e (0, 1) =^ 1 = A TT^^ + (1 - A) ttI with tt] , tt^ > tt] = = 1 . 

A similar argument can be applied for the entries of vr that are 0. It follows that we must have tt^ = tt = tt^, 
proving that tt is extremal. 

It remains thus to prove only that V has no other extremal points. To achieve this, it suffices to prove that 
the extremal points of V are all integer vectors, i.e. vectors all entries of which are integers - once this is 
established, the Lemma is proved, since the only integer vectors in V are those with all entries in {0, 1}. 

To prove that the extremal points of V are all integer vectors, we invoke the Hoffman-Kruskal theorem (see 
|15j . Theorem 7C.1), which states that, given a L x K matrix M, with all entries in {—1, 0, 1}, and a vector 
b e with integer entries, the vertices of the polytope defined by {x G M.^ ; {M.x)i <hi for ^ = 1, . . . , L} 
are all integer vectors in if and only if the matrix M is totally unimodular, i.e. if and only if every square 
submatrix of M has determinant 1, or —1. 



We first note that (C.l ) can indeed be written in this special form. The equality ^ Wij — M can be recast 
as the two inequalities j nij < AI and — j nij < ~M . The full system (C.l I can then be written as 



{M-k)i < bi tor £ ^ 1, ... ,L, where M is a (2iV + 2 + N"^) x matrix constructed as follows. Its first 2N 
rows correspond to the constraints on the sums over rows and columns; the entries of the next row are all 
1, and of the row after that, all —1 - these two rows correspond to the constraint X^i j "^ij — -^^5 ^he final 
jV^ X iV^ block is diagonal, with all its diagonal entries equal to —1. The first 2N entries of b are 1; the 
next 2 entries are M and — M; its final iV^ entries are 0. By the Hoffman-Kruskal theorem it suffices thus 
to show that M is totally unimodular. 

Because the last N'^ rows, the bottom rows of M, have only one non-zero entry, which equals —1, we can 
disregard them. Indeed, if we take a square submatrix of M that includes (part of) one of these bottom rows, 
then the determinant of the submatrix is if only zero entries of the bottom row ended up in the submatrix; 
if the one -1 entry of the bottom row is an entry in the submatrix, then the determinant is, possibly up to 
a sign change, the same as if that row and the column of the —1 entry are removed. By this argument, we 
can remove all the rows of the submatrix partaking of the bottom rows of M. 

We thus have to check unimodularity only for M', the submatrix of M given by its first 2N + 2 rows. If 
any submatrix contains (parts of) both the {2N + l)st and the {2N + 2)nd row, then the determinant is 
automatically zero, since the second of these two rows equals the first one, multiplied by -1. This reduces 
the problem to checking that M", the submatrix of M given by its first 2N + 1 rows, is totally unimodular. 

We now examine the top 2N rows of M" more closely. A little scrutiny reveals that it is, in fact, the adjacency 
matrix G of the complete bipartite graph with iV vertices in each partj^ It is well-known (see e.g. Theorem 
8.3 in [28]) that this adjacency matrix is totally unimodular, so any square submatrix of M" that does not 
involve the (2A^ -I- l)st row of M" is already known to have determinant 0, 1 or —1. Wc thus have to check 
only submatrices that involve the last row, i.e. matrices that consist of a (71 — 1) x n submatrix of G, with 
an added nth row with all entries equal to 1. We'll denote such submatrices by G'. 

We can then use a simple induction argument on n to finish the proof. The case n = 2 is trivial. In proving 
the induction step for n = m, we can assume that each of the top to — 1 rows of our m x m submatrix G' 
contains at least two entries equal to 1, since otherwise the determinant of G' would automatically be 0, 1 
or -1 by induction. 



The adjacency matrix A for a graph Q has as many columns as G has edges, and as many rows as Q has vertices; if we 
label the rows and columns of A accordingly, then A^e = 1 if the vertex v is an end point of the edge e; otherwise Aye = 0. An 
adjacency matrix thus has exactly two nonzero entries (both equal to 1) in each column. The number of nonzero entries in the 
row with index v is the degree of v in the graph. 
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The first to — 1 rows of C correspond to vertices in the bipartite graph, and can thus be partitioned into 
two sets Si and S2, based on which of the two parts of A'' vertices in the graph they pertain to. Let us call 
S the larger of 5*1 and 6*2; S consists of at least [^^^^1 rows. Let us examine the (#5') x to sub-matrix G" 
constructed from exactly these rows. We know that each column of G" has exactly one entry 1, since all the 
rows of G" correspond to the same group of vertices in the bipartite graph. Therefore, summing all the rows 
of G" gives a vector v of only zeros and ones; since each row in G" contains at least two entries equal to 1, 
the sum of all entries in v is at least 2 ([^^^^1) > to — 1. The vector v has thus at least to — 1 entries equal 
to 1; the remaining TOth entry of this linear combination of the top to — 1 rows of G' is either 1 or 0. In the 
first case, the determinant of G' vanishes, since its last row also consists of only ones. In the second case, 
we can subtract v from the last row of G' without changing the value of the determinant; the resulting last 
row has all entries but one equal to 0, with a remaining entry equal to 1. The determinant is then given by 
the minor of this remaining entry, and is thus 0, 1 or -1 by the unimodularity of G. □ 
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